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Abstract 

Using computer simulations and a thermodynamically self consistent integral equation we inves- 
tigate the phase behaviour and thermodynamic anomalies of a fluid composed of spherical particles 
interacting via a two-scale ramp potential (a hard core plus a repulsive and an attractive ramp) and 
the corresponding purely repulsive model. Both simulation and integral equation results predict 
a liquid-liquid de-mixing when attractive forces are present, in addition to a gas-liquid transition. 
Furthermore, a fluid-solid transition emerges in the neighbourhood of the liquid-liquid transition 
region, leading to a phase diagram with a somewhat complicated topology. This solidification at 
moderate densities is also present in the repulsive ramp fluid, thus preventing fluid-fluid separation. 

PACS numbers: 64.70.-p,61.20.Ja,05.20.Jj 



1 



I. INTRODUCTION 



The existence of a liquid-liquid (LL) equilibrium and density, diffusivity and structural 
anomalies in single component fluids has attracted considerable interest in the last decade. 
The density and diffusion anomalies present in liquid water (i.e. the existence of a maximum 
in the density within the liquid phase, and an increase in diffusion upon compression, up to 
a certain point) have long since been known, and have been reproduced by several of the 
existing interaction models^ 2 -^ 1 ^. Experimental evidence for LL coexistence has been found 
for phosphorous 6 , triphenyl phosphite^, and n-butanol^, and it has been suggested that this 
LL equilibrium might be the source of the anomalies encountered in water. Computer sim- 
ulations predict the existence of LL equilibria not only in water-^ but also in other loosely 
coordinated fluids, such as silicon 9 , carbon^ and silica 11 . Whether these transitions physi- 
cally exist or not is still open to debate, since in most cases they correspond to supercooled 
states which are rendered experimentally inaccessible by crystallisation. However, closely 
related first order transitions between low and high density amorphous phases have indeed 
been found for water—, silica^ and germanium oxide^. 

While the use of realistic models can provide reasonable explanations for the experimental 
behaviour of a physical system, a more thorough description of the mechanisms underlying 
both LL phase transitions and density, diffusivity and related anomalies can be acquired 
via the study of simplified models. In the case of LL equilibrium and polyamorphism in 
molecular fluids^ the simple model of Roberts and Debenedetti^ has successfully accounted 
for the behaviour of network forming fluids^. On the other hand, since the pioneering work 
of Hemmer and Stell^^ 9 - it is known that a simple spherically symmetric potential in which 
the repulsive interaction has been softened (in this case by the addition of a repulsive ramp 
to the hard core) can lead to the existence of a second (LL) critical point, as long as a 
first (liquid-vapour, LV) critical point existed due to the presence of a long range attractive 
component in the interaction potential. As well as the ramp potential, other simple models 
with two distinct ranges of interaction, such as the hard-sphere square shoulder- square well 
potential studied by Skibinsky et al.— , have also been shown to exhibit LL equilibria. Indeed 
the presence of two interaction ranges explains the competition between two locally preferred 
structures (LPS) - a LPS being defined as an arrangement of particles which, for a given 
state point, minimises some local Helmholtz energy 21 . This competition between two LPS 
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helps to rationalise the existence of polyamorphism and LL equilibria in single component 
glassy systems and fluids^. 

More recently, the original model proposed by Hemmer and Stell has regained attention, 
especially since Jagla3^ stressed the similarities between the behaviour of the Hemmer-Stell 
ramp potential and the anomalous properties of liquid water. This has been further explored 
by Xu et al.- 1 ^ who analysed the relationship between the LL transition and changes in 
the dynamic behaviour of fluids interacting via a soft core ramp potential with attractive 
dispersive interactions added. Gibson and Wilding- 4 have recently presented an exhaustive 
study of a series of ramp potentials exhibiting LL transitions and density anomalies, whose 
relative position and stability with respect to freezing might be tuned by judicious changes 
in the interaction parameter. Moreover, Caballero and Puertas^ have also focused on the 
relation between the density anomaly and the LL transition for this model system by means 
of a Monte Carlo based perturbation approach. The aforementioned authors find that, in 
this case, the density anomaly is absent when the range of the attractive interaction is 
sufficiently small. 

In this work we shall refer to this interaction potential as the 'attractive two-scale ramp 
potential' (A2SRP). Additionally, when one considers the system without any attractive 
contribution, i.e. a hard-sphere core plus a repulsive ramp, as the 'two-scale ramp potential' 
(2SRP). Although it has been found that the system exhibits density anomalies, it is a pos- 
sibility that LL transition is preempted by crystallisation^. In this purely repulsive model, 
the relation between static and dynamic anomalies has been explored in detaiP^- 1 ^ (simi- 
lar results have been found for a dumbbell fluid with repulsive ramp site-site interactions^) 
and the connection between these anomalies and structural order has been investigated by 
Yan et al.— with the aid of the Errington-Debenedetti order map^ 1 ^. These authors found 
that there is a region of structural anomalies (in which both translational and orientational 
order decrease as density is increased) that encapsulates the diffusivity and density anomaly 
region. The intimate relationship between transport coefficients, and the excess entropy was 
made clear first by Rosenfeld^ 1 ^, and again later, specifically for the atomic diffusion, by 
Dzugutov^. Needless to say, a corollary of this relation is that any anomalous diffusion will 
be accompanied by an anomalous excess entropy. Recently this link has been shown to hold 
for both liquid silica and for the 2SRP model by Sharma et al.— , and for the discontinuous 
core-softened model by Errington et al.— . 
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The principal objective of this work is to extend our knowledge of the phase behaviour 
of systems interacting via either attractive or purely repulsive two scale ramp potentials 
(A2SRP and 2SRP). To that purpose exhaustive Monte Carlo calculations have been per- 
formed in order to determine the phase boundaries of the gas, liquid and solid phases for 
both model systems up to moderate densities -slightly beyond the high density branch of 
the LL equilibrium. Self-consistent integral equation calculations performed on the A2SRP 
model complement the Monte Carlo results and agree qualitatively as to the location of the 
LL equilibria and quantitatively for the LV equilibria. The location of Widom's line (the 
locus of heat capacity maxima) is obtained for both models. The temperature of maximum 
density (TMD) curve is obtained for the repulsive 2SRP model and correlated with the 
location of Widom's line. The rich variety of phases present in these simple models will be 
illustrated in the calculated phase diagrams. 

The structure of the paper is as follows; the model and computational procedures used 
herein are introduced in Section II. In Section III the most significant results are presented 
and discussed. Finally, the main conclusions and future prospects can be found in Section 
IV. 



II. MODEL AND COMPUTATIONAL METHODS 



The first model system consists of hard spheres of diameter a, with a repulsive soft core 
and an attractive region. The interaction potential reads 



u{r) 



oo if r < o 

W r -{W r -W a )j^ ifa<r<d a 

if d a < r < d c 
if r > d r 



(i) 



o 



where W r > and W a < 0. The values used here for the parameterisation of the model 
are the same as those used in Ref. []. Thus, energy units are defined by \W a \, with the 
reduced temperature T* = k B T/\W a \ {ks being the Boltzmann constant) and we have set 
W r /\W a \ = 3.5. The units of length are reduced with respect to the hard core diameter, and 
so one has d a /a = 1.72, d c /a = 3.0 and the reduced density is p* = per 3 as usual. This set of 
parameter values and Eq.(TjQ) define the A2SRP model. For the purely repulsive system, we 
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FIG. 1: Attractive two scale ramp (A2SRP, solid line) and repulsive two scale ramp (2SRP, dash- 
dotted line) potential models. 



have chosen the corresponding repulsive potential that would result from a Weeks-Chandler- 
Anderson^ decomposition of Eq.([T]), namely 



u r (r) = < 



oo if r < a 

W r -W a - (W r - W a )^ iia<r<d a 
if r > d a 



(2) 



In Figure [T] both the attractive and purely repulsive interactions are illustrated. 



A. The self consistent integral equation approach 

The integral equation calculations are based on the Ornstein-Zernike relation, which for 
simple fluids reads 3 -^ 

h(r 12 ) = c(r 12 ) + p J c(r 13 )h(r 32 )dv 3 (3) 

where h(r) is the total correlation function (related to the pair distribution function, g, by 
h = g — 1) and c(r) is the direct correlation function. This equation requires a supplementary 
relation, whose general form is 38 

h(r) = exp [—/3u(r) + h(r) - c(r) + B(r)) - 1 (4) 

with (3 = {ksT)~ l . In Eq.Q the bridge function, B(r), is a diagrammatic sum of convolu- 
tions of h(r), and must be approximated. The simplest instance is the hyper-netted chain 



approximation (HNC), which implies B(r) = 0. Interestingly, this approximation predicts 
the existence of a LL equilibrium in the square shoulder- square well model studied in Ref. 139 
whereas for the A2SRP model only the LV equilibrium is reproduced. Moreover, in the case 
of the 2SRP model, it completely misses the density anomaly 26 . Preliminary calculations 
with a more elaborate closure such as the one proposed by Martynov, Sarkisov and Vompe 
(MSV)^, show that although it turns out to be extremely accurate in the determination 
of the LV equilibria, it fails to capture the LL transition. However, one observes a second 
Widom line at moderate densities, which seems to indicate an anomaly in the region where 
the LL is expected to appear. This implies that thermodynamic consistency should play a 
central role if the LL transition or the density anomaly is to be found. This is confirmed by 
the results Kumar et al.— for the 2SRP model. A fairly successful self consistent approach 
is the so-called Hybrid Mean Spherical approximation (HMSA) which smoothly interpo- 
lates between the HNC and the Mean Spherical Approximation (MSA) closures^i^. The 
corresponding closure reads 

g(r) = exp(-/?w, r (r)) 

v A , exp[/(r)Q(r) - c(r) - (3u a {r))} - l \ 

y \ l + W) J (5) 

with the interpolating function f(r) = 1 — exp(— ar). The repulsive component of the 
interaction given by Eq.([2]), and the attractive component simply given by u a (r) = u(r) — 
u r (r) with u(r) having been defined in Eq.(TjQ). Note that for the purely repulsive system, 
where u a = 0, one recovers the Rogers- Young (RY) closure^. The parameter a is fixed 
by requiring consistency between virial and fluctuation theorem compressibilities. With the 
virial pressure defined by 

(3P V = P ~ f P 2 T r*^ g (r)dr + f P V^+) (6) 



and 



(3k t x /p 



df5P 



dp 

Consistency could be implemented on a local level by equating 

d(3P v 



= 1 — Airp / r c(r)dr. (7) 

T JO 



(3k t 1 /p 



(8) 



dp 

Alternatively, one might use a global consistency. In this case, in the absence of spinodal 
lines, one can simply enforce the condition 

(3P v (p,T) = pP K (p,T) (9) 
6 



with 



Anp / r 2 c(r; p', T)dr 



dp'. 



(10) 



In the presence of spinodal lines one might use a mixed thermodynamic integration path, as 
suggested by Caccamo et al.— . In this paper extensive use is made of the local consistency 
(LC) approach. We have also explored the use of global consistently (GC) as an alternative 
to overcome the deficiencies exhibited by the LC approach. From a practical point of view, 
in order to calculate the derivative of the virial pressure in Eq.([8]), we shall assume that the 
parameter a is locally independent of the density, as is customary. For density independent 
potentials this leads to minor discrepancies when comparing the integrated pressures in Eq. 
([H]), so both LC and GC approaches yield similar results 4 ^. Here the situation is somewhat 
different, as will be seen in the proceeding Section. Finally, the right hand side of Eq.([H]) is 
explicitly calculated via 



d(3P v 



dp 



- 1 



2(3P V 
P 



2vr 2 3 d 9(r) 
3 P 



du r {r) \ ( dg{r) 



dr 



dp 



dr 



dp 



(11) 



where the density derivative of the pair distribution function is obtained by iterative solution 
of the integral equation that results from differentiation of ([3]) and ([5]) with respect to 
density. This procedure, first proposed by Belloni 4 ^ is substantially more efficient than the 
classical finite differences approach when it comes to evaluating the derivative of f3P v . In 
addition to the density derivatives of the pair correlation function we have also evaluated 
the corresponding temperature derivatives in order to calculate the heat capacity 47 



C* 



Nk 



B 



3 
2 



2vrp / (3u(r)g(r) f3u{r) - f3 



duj{r) 



8(3 



r dr 



(12) 



where u = h — c + B. Note that the integral equations in terms of the derivatives of the 
correlation functions (either with respect to p or (3) are extremely well conditioned and are 
rapidly solved, even when starting from ideal gas initial estimates. 
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B. Simulation algorithms 

Initially the principal aim was to establish whether the A2SRP model does indeed exhibit 
a fluid-fluid equilibrium. Therefore, in order to compute both the liquid-vapour and the LL 
equilibria, a procedure^ based on the algorithm of Wang and Landau^* 5 ^ was implemented. 
Subsequently, it was necessary to compute the equilibrium of the fluid phases with the low 
density solid phase. To that end use was made of thermodynamic integration 51 and the so- 
called Gibbs-Duhem integration techniques 5 ^* 5 ^ 1 ^, which took advantage of the behaviour 
of the repulsive model in the low temperature limit. The use of two different methods in 
independent calculations provides a check on the quality of the results. 

1. Wang-Landau-type algorithm 

A more detailed explanation of the Wang-Landau procedure can be found elsewhere^. 
Here, we shall restrict ourselves to outlining its most salient features. The algorithm fixes the 
volume and temperature and then samples over a range of densities (N) of the system. For 
each distinct density the Helmholtz energy function F(N\V, T) is calculated. The procedure 
resembles, to some extent, well known grand canonical Monte Carlo methods 5 ^* 5 ^, but using 
a different (previously estimated) weighting function for the different values of N. Once 
F(N\V, T) is known for different values of N, one is able to compute the density distribution 
function for different values of the chemical potential, /x, and extract the conditions for which 
the density fluctuations are maximal. This could indicate the presence of phase separation 
and eventually, by means of finite-size-scaling analysis, one can locate regions of two phase 
equilibrium and the locus of critical points.— By following such a procedure one can easily 
compute the liquid-vapour (LV) equilibrium at high temperatures. The computation of the 
liquid-vapour equilibrium was predominantly performed with a simulation box of volume 
V/a 3 = 1000 (i.e. box side length of L* = L/a = 10). In order to obtain a precise location 
of the critical point simulation runs were performed close to the estimated critical point for 
different system sizes: L* = 6,7,8,9,10,11, and 12. By means of re-weighting techniques 
and finite-size-scaling analysis, we have estimated the location of the critical point to be: 
T* = 1.487 ± 0.003, p c a 3 = 0.103 ± 0.001, and p c o- 3 /\W a \ = p* ~ 0.042. 

At low temperatures the LL equilibrium was found, but as the temperature was further 
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reduced we were unable to obtain convergence in the estimates of the Helmholtz energy 
function, and the systems eventually formed a low density solid phase. Given this situation 
it was surmised that a triple point, if it exists, between the vapour and the two liquid phases, 
could not be stable, and the presence of a stable low density solid phase is most likely. After 
performing simulations at several temperatures, with different system sizes (0 < pa 3 < 0.55, 
L* = 6, 8, 10), it was found that the LL equilibrium appeared for temperatures below T* ~ 
0.38. Using the results for L* = 10, we can estimate the critical point for the LL equilibrium 
at T* ~ 0.378 ± 0.003, p c a 3 ~ 0.380 ± 0.002 and (p* c /T*) ~ 0.49 ± 0.01. 

2. Thermodynamic Integration: 2SRP model 

In order to compute the phase equilibria between the low density solid and the fluid 
phases, it is useful to firstly consider the repulsive model. Jagla3^ computed the stability of 
different crystalline structures at T = for the repulsive model supplemented by a mean field 
attraction term, for a ramp model with d a ja = 1.75, and found that the stable crystalline 
phases at low pressure correspond to the face-centred cubic (fee) and hexagonal close packed 
(hep) structures. On the other hand, in the limit T — > 0, and not too high a pressure the 
2SRP model approaches an effective hard sphere system with diameter d a . We can, then, 
expect the fee structure to be the stable phase in the ranged : 1. 038 (a/d a ) 3 < pa 3 < 
\^2(a/d a ) 3 , with the transition from the low density fluid to the solid at a pressure^ of 
j3pa 3 ~ 11.6 (a/d a ) 3 . Assuming that the fee is indeed the structure of the low density solid 
phase, and taking into account the low temperature limit for the 2SRP model one can use 
standard methods to study the equilibrium of this solid with the fluid phase(s). 

In order to relate the two ramp models considered in this work, and to explain the 
procedures used to compute the phase diagram of the systems using thermodynamic and 
Gibbs-Duhem integration, it is useful to write down a generalised interaction potential, Ui(r) 
using an additional 'perturbation' parameter, A: 

Ui(r\X) = u r (r) + \u a (r). (13) 

When A = 0, there are two interesting limits; when T — > oo the model approaches that of a 
system of hard spheres of diameter a, whereas in the limit T —* one once again encounters 
a hard sphere system, but now having a diameter of d a . In order to compute the Helmholtz 
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energy for condensed phases of the ramp potential systems, one can perform thermodynamic 
integration over the variable (3 and the parameter A. Thus one is able to compute the 
difference between the Helmholtz energy of a particular state, and the aforementioned hard 
sphere systems. 

In the canonical ensemble the partition function can be written as: 

yN r 

Qi = J dR*exp[-/3EMV,R*,A)] (14) 

where A, the de Broglie thermal wavelength, depends on the temperature; R* represents 
the reduced coordinates of the N particles: (J dR* = 1), and U\ is the potential energy 
(given by the sum of pair interactions Ui(r\X)). The Helmholtz energy function is given by: 
(5Fi = — logQi- One can then write the Helmholtz energy per particle, / = F/N, as: 

Pf = Pf id + Pf ex ; (15) 

where the ideal contribution f ld can be written as 

/3f d = 31og^ + log(pa 3 )-l. (16) 

<7 



and the excess part as 



Hr = ~ log { I dR* exp [-PUM R*, A)] }. (17) 



N 

In what follows we shall consider the Helmholtz energy as a function of the variables N, V, T, 
and A (or, equivalently, as N, p, (3 and A). We can then compute the derivatives of f3f ex with 
respect to (5 and A. I.e. 

and 

where the averaged quantity, (U a ), is given by: 

/N-l n \ 

w = EE ) ■ ( 2 °) 
\ i=\ j=i+i I 

On the other hand, using the well known thermodynamic relation: 



dF 
dV 



= ~P (21) 

NTX 
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one obtains 

d(Pf) _ PP 

dp ~ f {ZZ) 
For the 2SRP model (A = 0), one can compute the Helmholtz energy of the low density 
solid using thermodynamic integration. As mentioned earlier, in the limit T — > one has an 
effective hard sphere system of diameter d a . The phase diagram of the hard sphere system^!, 
the equation of state for the fluicr^ 1 ^ and the equation of state for the solid^ are all well 
known. Thus 

Pf ex (p, T) = m:(pdl) - £ dT x . (23) 

where /3fHg(pd 3 ) is the excess Helmholtz energy per particle in the /cc-solid phase and U/N 
is the excess internal energy per particle. The integrand in (1231) is well behaved in the limit 
T — > 0. Simulations have been undertaken for iV = 500, pa 3 = 0.24 (pd 3 ~ 1.2212), for 
various temperatures; T* = 0.025, 0.050, 0.10, 0.15, 0.20, 0.45, 0.500. At T* = 0.500 
the system melts. The results for U /T 2 as a function of T have been fitted using a polynomial. 
Thus we are able to obtain a function to compute the Helmholtz energy of the low density 
solid phase in the range < T* < 0.45 at the reference reduced density of p^<J 3 = 0.24. 

In order to acquire data for the calculation of the fluid-low density solid equilibrium 
a number of simulations were performed for the low density solid phase, along several 
isotherms, having densities around p$ \ The low density solid is stable only within a small 
region of the p — T phase diagram, as can be seen in Fig. 10. The pressures in this stability 
range can be fitted to a polynomial, 

Pp s (T,p) = J2* ( i\T)p\ (24) 

i=0 

We can then compute the Helmholtz energy of the low density solid, as a function of the 
density on each isotherm, using 

(3f( Pl , T) = pf(po, T) + f 1 ^^-dp, (25) 

where p = p { s \ and /3p(T,p) = (3p s (T,p). 

The Helmholtz energy of the low-density fluid phase is computed using the results of 
several simulations, having = 500, at several temperatures and densities; typically pa 3 = 
i x 0.025, with i = 1,2..., 10. In order to guarantee that the samples were indeed within the 
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fluid phase the simulation runs were initiated from equilibrated configurations of temperature 
T* ~ 2. The pressure was fitted to a virial equation of state: 
PPg(P,T) ~ p 



P 2 



B 2 (T) + B 3 (T)p + B 4 (T)p 2 + ■■■, (26) 



which is then used in the calculation of the excess Helmholtz energy of the gas: 

t>r(P,T)=[' ' Mpi -V-» d» (27) 

Jo Pi 

Using (I15|I16II26||2~TI) one can compute the Helmholtz energy of the gas phase as a func- 
tion of the density. The Helmholtz energies of the high density fluid can be computed 
using thermodynamic integration from the high temperature limit. Performing a number of 
simulations, with N = 500, for different values of (3 at a fixed reference density, pi, we have 

Pf ex (pi, 0) = /3/LW) + J P U(N ^ Pl) dp, (28) 

For selected temperatures several simulation runs are performed for various densities, and 
the results then are used to fit the equation of state of the liquid branch, 

(3p{T,p) = Y,a i t\T)p\ (29) 

i=0 

Similarly, we can compute the Helmholtz energy of the high density fluid at a given tem- 
perature by again performing simulations at several densities, and then fitting the results of 
the pressure as a function of the density and using ( 1251) with p = /3q where p[, cr 3 = 0.40. 
Once we have the equations of state for the low density fluid, the low density solid, and the 
high density fluid, and the corresponding Helmholtz energies at a certain reference density 
for each case, it is straightforward to compute the equilibria by locating the conditions at 
which both the chemical potential and pressure of two phases are equal. 



3. Thermodynamic Integration: A2SRP Model 

The Helmholtz energy of the fluid phases can be computed following the same procedures 
outlined in the previous subsection. In order to compute the Helmholtz energy of the low 
density solid phase one can perform thermodynamic integration starting from the repulsive 
model and integrating f|T9|) at constant /3, p, and N: 

(3 f a (p, T) = (3f r (p, T) + ^ j\u a (N, p, T, A)) dX, (30) 
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where the indices a and r indicate the A2SRP and the 2SRP models respectively. Such 
an integration was carried out using A = k x 0.10 (with k — 0, 1, 2, ■ • • , 10) with N = 500 
at pa 3 = 0.24 at the temperatures T* = 0.30 and T* = 0.35. The integrand (U a ) is well 
behaved in both cases. The Helmholtz energy of the A2SRP model low density solid phase 
at the reference density p^cr 3 = 0.24 as a function of T was parametrised via simulation 
results. Integration of (1301 at two distinct temperatures provided a check of the numerical 
consistency. As for the 2SRP model various simulations were performed in the region of p^ 
in order to compute the equation of state and the Helmholtz energy of the low density solid. 

Thermodynamic integration was used to study the equilibrium between the low density 
solid and the high density liquid. The procedure was analogous to that used to study the 
low density solid-high density fluid equilibrium of the 2SRP model . To check the results of 
the LL equilibrium thermodynamic integration was performed at T* = 0.35. In order to do 
this the Helmholtz energy was calculated at the 'reference' reduced densities of 0.30 and 0.45 
using the procedure derived from (T281) . and then by performing a simulation at T* = 0.35 
at several densities. The results for the equation of state were fitted for both branches, and 
the conditions of thermodynamic equilibria were calculated. A good agreement (within the 
error bars) was found with the results from the Wang-Landau calculation. 



4- Gibbs-Duhem Integration: 2SRP Model 

The partition function, Q(N,p,T, A) in the isothermal-isobaric (NpT) ensemble can be 
written as^i: 

Qn p t = /3p j dV exp [-(3pV\ Qnvt(N } V, T, A). (31) 
The chemical potential, p is given by 

(3p = -^\ogQ NpT . (32) 

Considering constant N, and p = p(p, T, A) one has: 

d((3p) = ^d{3 + ¥jd(J3p) + ^d\. (33) 

where E is the internal energy (including kinetic and potential contributions). Let us con- 
sider two phases, a and (3, in conditions of thermodynamic equilibrium (i.e. equal values 
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of T, p and /i); if one makes a differential change in some of the variables j3, (ftp) and A, 
thermodynamic equilibrium will be preserved if: 



d[A(Pfj)] = = AUdf3 + AVd(pp) + AU a d\. (34) 

where AX expresses the difference between the values of the property X in both phases, 
and X = X/N. In (134]) we have considered equal values of the kinetic energy per particle 
in both phases (classical statistics). 

In order to compute the fluid-low density solid equilibrium of the 2SRP model we make 
use of the Gibbs-Duhem integration scheme. As the starting point of the integration we 
consider the fluid-solid equilibrium of the effective hard sphere system, with diameter d a at 
T = 0. We have integrated (134")) for A = 0. After a number of short exploratory simulation 
runs, the Gibbs-Duhem integration computation for the 2SRP model was performed in three 
subsequent stages: i) low density fluid-low density solid equilibrium at low temperature, ii) 
fluid-low density solid at 'higher' temperatures, and iii) low density solid-high density fluid 
at low temperatures. For the first stage equation (1511) was modified to obtain (in reduced 
units) the finite interval numerical approach: 

<»> 

The calculation of U and V for both phases at the estimated coexistence conditions was 
performed by Monte Carlo simulations with N = 500, in the NpT ensemble. At T* = 
the coexistence of the fluid and the solid phases occurs 51 at (p*/T*) ~ 2.286 (i.e. ~ 11.6 x 
(a/d a f). A number of simulations (T* = 0.01,0.02 with (p*/T*) ~ 2.29) were used to 
estimate the initial values of the slope d(p*/T*)/dT*. The integration was then carried out 
from T* = 0.01 to T* = 0.40 with a step of 5T* = 0.01, 

In the second stage of the integration the coexistence temperature reaches a maximum, 
thus the independent variable in the integration was substituted for 

T 2 AV 

ST^^^dip/T). (36) 

This integration was started from the pressure at which the temperature reaches T* = 
0.40 (i.e. (p*/T*) ~ 3.65) using an integration step of 5(p*/T*) = 0.05. Initially, the 
temperature increases with p until the coexistence line reaches a maximum temperature; at 
this point the density of both phases is equal: (T = 0.422 ± 0.002, p* /T* = 4.65 ± 0.05, 
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per 3 = 0.260 ± 0.002). For higher pressures the fluid density becomes higher than the low 
density solid density, and the temperature decreases with p/T. The third stage is initiated 
upon reaching T* = 0.40; which happens at (p*/T*) ~ 6.33. This involves reverting to the 
integration scheme given by (|35p . this time using 5T* = —0.01. 

5. Gibbs-Duhem Integration: A2SRP Model 

In order to compute the fluid-low density solid equilibria of the A2SRP system, this 
system is connected to the fluid-low density solid equilibria of the 2SRP model, by tuning 
the perturbation parameter A from A = to A = 1, while maintaining the conditions 
of thermodynamic equilibrium, i.e. those of Eq. f|34|) . For this system this connection 
stems from the initial state: A = 0, T* = 0.34, (p*/T*) = 3.05, which corresponds to the 
equilibrium between the low density solid with a lower density fluid. This was done by 
keeping T constant and integrating numerically the equilibrium pressure as a function of A, 
i.e.: 

5(p*/T*) ~ (37) 

using SX = 0.025. The final point was found to be (A = l,p*/T* = —0.324), corresponding 
to a equilibrium between a low density solid with density p* = 0.248 and a low density liquid 
with p* = 0.239. Given the negative sign of the pressure the equilibrium found is metastable. 
From this point we can retake the Gibbs-Duhem integration (now for the attractive model, 
i.e. setting A = 1); using fl36l) . with an integration step of 5(p* /T*) = 0.01. As with the 
repulsive model a temperature maximum was found at Tfa ~ 0.345, per 3 ~ 0.255. However, 
in this case such a maximum seems to be metastable: (p*/T*) ~ —0.01. The solid-fluid 
equilibrium becomes thermodynamically stable at p*/T* ph 10 -4 (where this coexistence 
line crosses the liquid- vapour equilibrium line). At higher pressures and lower temperatures 
this low density liquid-low density solid coexistence line will eventually coincide with the LL 
equilibrium line at a triple point. Beyond this a stable low density solid-high density liquid 
equilibrium will appear. 

In order to obtain a precise estimate of the conditions for which the two liquid phases and 
the low density solid are in equilibria (i.e. the location of the triple point), Gibbs-Duhem 
integration was performed for the LL equilibrium. Both thermodynamic integration and the 
Wang-Landau procedure were used to calculate the LL equilibrium, at T* = 0.35. Both of 



15 



the results agreed to within the error bars: (j3pa 3 ~ 0.465, with reduced densities of the 
liquid phase: pier 3 ~ 0.30, and P2O" 3 — 0.46). From this initial point the Gibbs-Duhem 
integrations of the phase equilibrium were carried out using (155]) with ST* = 0.0025. The 
LL equilibrium line met the low density liquid-low density solid equilibrium at T* ~ 0.331, 
(p*/T*) ~ 0.437, with the three phases having reduced densities: p(i ow density solid) ' 3 — 0.263, 
P(iow density iiq.)°" 3 — 0.293 and P(hi g h density Hq.) " 3 — 0.48 . Gibbs-Duhem integration was then 
used to calculate the low density solid-high density liquid equilibrium, starting from this 
triple point, using (|35|) with 6T* = 0.01 

6. Vapour-solid equilibrium 

The computation of the density of the low density solid at equilibrium with the vapour 
phase for the A2SRP model is straightforward. In practice, due to the low pressure at which 
this equilibrium occurs, the NpT simulations of the low density solid are undertaken at 
ftp ~ 0, again with iV = 500. This immediately yields an accurate estimate of the density 
of the solid in equilibrium with the vapour phase. 

7. Details of the Gibbs-Duhem integration scheme 

A simple predictor-corrector scheme was used to build up the coexistence lines. Generi- 
cally one seeks the solution of: 

£ = /(*,»). (38) 
from some initial condition x , y . In the present case the function f(x, y) has to be computed 
via a pair of computer simulations. This result will be affected by statistical error. Using a 
set of discrete values Xi = xo + i ■ h, the following scheme was used: 

yffi = Vi + I m - /i-i) h (39) 

Vi+l =Vi J r\{fi J r fi+l) ■ (40) 

with: 

fi = f(x i ,y? ) ) (41) 

where ( 139]) and (1401) are respectively the predictor and corrector steps. This simple algorithm 
is both accurate and robust. 
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1 2 3 4 5 

r/a 

FIG. 2: (Colour online) Pair distribution function for the A2SRP fluid at low (pa 3 = 0.10) and 
high (pa 3 = 0.80) densities at T* = 1.8. Inset: magnified view of the region r/a from 2.5 to 6. 

III. RESULTS 

In order to illustrate the two competing LPS that the A2SRP model system exhibits at 
low and high densities (or low and high pressures respectively) in Figure [2] we have presented 
the pair distribution functions for two representative states of the A2SRP model, obtained 
from the HMSA integral equation. The high density state is close to a hard sphere fluid of 
diameter a, whereas the low density g(r) corresponds to a fluid of soft particles of diameter 
d a = 1.72a. The purely repulsive system leads to similar results, with only minor differences, 
due to the lack of dispersive forces (the sharp minima at 3a are absent). It is clear that the 
LL transition will result from the 'chemical equilibrium' between two essentially different 
'fluids' that stem from the same interaction. 

A. Liquid-vapour and liquid-liquid equilibria 

In order to calculate the phase diagram using the HMSA integral equation, one has 
to resort to thermodynamic integration. Although direct closed formulae for the chemical 
potential are available in the work of Belloni^, it has been the authors experience that better 
results are obtained when using thermodynamic integration. For a given state point, at 
density p (on the right hand side of the phase diagram) and a subcritical inverse temperature 
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FIG. 3: (Colour online) Phase diagram of the A2SRP model as obtained from the HMSA and 
Wang-Landau Monte Carlo simulations. HMSA results for the loci of maxima of C* (the Widom 
line), isothermal compressibility (j>kt/P) and the set of thermodynamic states for which the first 
non trivial minimum of h(r) vanishes. Additionally, the boundaries of the non-solution region of 
the integral equation and the thermodynamic spinodal of the LL transition are also depicted. 




FIG. 4: (Colour online) Inverse virial and fluctuation theorem compressibilities (solid and dashed 
curves respectively) as calculated using the HMSA integral equation at T* = 0.38 for the attractive 
model. Note how the consistency degrades when approaching the LL transition region (see Figj3]) 
to the point that the virial compressibility crosses a thermodynamic LL spinodal before hitting the 
LV spinodal, whereas the fluctuation theorem compressibility only detects the presence of the LV 
spinodal. 
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(3, one can obtain the free energy per particle along a mixed path of the form 

fW) =iJ p mp',Po)-p)/p /2 d P > 



Po 

+i J* U(p, (3')/Nd(3' + log P A 3 - 1 (42) 



where (3q is a supercritical inverse temperature With this expression one is able to evaluate 
the chemical potential by using the relation 

P 

The phase equilibria can be determined by equating pressures and chemical potentials for 
both the gas and liquid phases, and for both the low density liquid and high density liquid 
phases. 

The phase diagram thus obtained is plotted in Figure [3] together with the corresponding 
Monte Carlo estimates. Additionally, the curves corresponding to the loci of maxima of 
C* (the Widom line), maxima of the isothermal compressibility (pKy//?) and the set of 
thermodynamic states for which the first non trivial minimum of h(r) vanishes are presented. 
This latter quantity separates supercritical states with gas-like local order from those with 
liquid-like order- 8 . Interestingly, as found in the Lennard- Jones system with a completely 
different closure 58 , these three singular lines are seen to converge towards the LV critical 
point. Although this should be expected from the locus of heat capacity and compressibility 
maxima, there is no apparent fundamental reason why this should also happen for the line 
of states for which the first minimum of h(r) = 0. Moreover, we have found that similar 
results (with different locations for the critical point) are obtained using different closures, 
such as the HNC, Reference HNC or the MSV. Finally, one sees that the upper part of the 
binodal line is broken, since convergent solutions cease to exist before the critical point is 
reached. In fact, a small portion of the high temperature LV binodal has been obtained 
by extrapolation of the chemical potentials and pressures at lower vapour densities (for 
this reason the no-solution and the binodal curves cross). This is a well known limitation 
of a number of integral equation approaches, encountered even when treating the simple 
Lennard- Jones fluid^. 

As for the LL transition, in Figure [3] the HMSA binodal results obtained from the ther- 
modynamic integration, its corresponding thermodynamic spinodal curve and the Wang- 
Landau Monte Carlo binodal estimates are presented. A second Widom line that crosses 
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FIG. 5: (Colour online) Residual inconsistency between virial and fluctuation theorem pressure as 
a function of the mixing parameter a in the HMSA closure (Eq. [5]) for the 2SRP model. 

the LL binodal in the neighbourhood of the HMSA critical point is also found. As well as 
this, a second line of isothermal compressibility maxima seems to indicate the presence of a 
LL transition, approaching the simulation LL critical point. Note, however, that in this case 
the no-solution line does not provide any clue as to the presence of a phase separation. In 
contrast to the HMSA behaviour in the vicinity of the LV critical point, now the curves of C* 
and pkt//3 maxima do not converge to a common estimate of the critical point. This reveals 
a substantial inconsistency between the virial and the fluctuation theorem thermodynamics. 
This is best illustrated in Figure H] where the virial and fluctuation theorem compressibili- 
ties are plotted. Here it is observed that the local consistency criterion of Eq.([HD leads to 
a good agreement for densities pa 3 > 0.4, but a considerable inconsistency shows up at 
lower densities, corresponding to the LL transition (see Figj3]). The assumption of a density 
independent a parameter in the HMSA closure (Eq. [5]) lies at the root of this failure. Thus, 
whereas the virial pressure predicts a LL transition, no diverging correlations appear near 
the thermodynamic spinodal, and the equation has solutions throughout the two phase re- 
gion. This is in marked contrast to the results for the LV transition. On the other hand, 
in FigJH the large values of (d(3P/ dp)T in the region 0.2 < pa 3 < 0.3 show that the fluid 
is almost incompressible, which is somewhat surprising, especially at these relatively low 
densities. This indicates that a solid phase (or a glassy state) may well be about to appear. 
We have attempted to explore alternatives to this breakdown in the LC approach. After 
an unsuccessful attempt to incorporate a linear density dependence on the a parameter in 
the HMSA interpolating function, in conjunction with a two parameter optimisation, we 
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FIG. 6: (Colour online) Zero pressure isobar calculated using the HMSA scheme and also by means 
of NPT MC simulation, for the A2SRP model. 

focused on the GC alternative expressed in Eq.([9]). This can be easily implemented using 
the strategy proposed by Caccamo et al. 44 . To make matters simpler, the GC HMSA calcu- 
lations were restricted to the purely repulsive 2SRP model. At relatively high temperatures 
(T* « 2) the procedure worked well and lead to results slightly better that those of the 
LC HMSA. However, as the temperature is lowered, one finds that the optimisation loop 
diverges even for relatively low densities. The reason behind this divergence is illustrated 
in FigJHl It was observed that for two neighbouring states, a slight decrease in temperature 
upwardly displaces the residual inconsistency curve. As a consequence of this, when looking 
for consistency, the minima are now found to be greater than zero, thus the numerical itera- 
tions will lead to either a divergence or result in an oscillating behaviour. This also explains 
why the alternative procedure based on the implementation of a density dependent a failed 
as well. A completely parallel situation is to be found for the attractive potential model. 
It is now clear that a more accurate treatment of this type of soft core models requires not 
only the implementation of GC conditions on the pressure, but also the use of more flexible 
closure relations. 

Finally, a few words regarding the comparison of the integral equation estimates with 
the simulation results. As far as the LV transition is concerned, the HMSA predictions are 
satisfactory. The re-entrant behaviour that is clearly seen on the high density side of HMSA 
LV curve is also found in the MC estimates, and can be more clearly appreciated in the zero 
pressure isobar depicted in Fig. El Here we see that the errors in the density do not exceed 
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FIG. 7: (Colour online) TMD curves as determined by HMSA calculations and NPT MC simulation 
for the purely repulsive model. Curves of maxima in the isothermal compressibility and heat 
capacity (Widom's line) calculated using the HMSA closure. 



5 per cent. This partly reentrant behaviour is a characteristic indication of the proximity 
of a triple point, and has also been found in various water models^i. On the other hand, 
whereas the LL critical temperature is reasonably well captured by the theory, the critical 
density is substantially underestimated, and the global shape of the LL equilibrium curve is 
not well reproduced. This can certainly be ascribed to the poor consistency of this HMSA 
approach in this region. Note, however, that this is the only integral equation theory, of 
those that we have checked, that captures the presence of the LL transition. 

When the attractive interactions are switched off, the LV transition disappears. The 
HMSA calculations yield lines of compressibility and heat capacity maxima that could be 
interpreted as indications of a possible LL transition (see Figure CJ). These lines appear in 
approximately the same locations as in the attractive model (see Figj3]), which indicates that 
this feature is entirely due to the presence of the soft repulsive ramp in the interaction. The 
possibility of a LL transition in this model has already been speculated upon by Kumar et 
al.— on the basis of the low temperature asymptotic behaviour of a series of isochores. The 
same conclusion might well be drawn from this work; however, once again calculations at 
lower temperatures are hindered by lack of convergence. It shall be seen in the next section 
that the fluid-solid equilibrium preempts the LL phase separation. 



As in Ref. |26|, we are able to calculate the TMD curves (lines for which (dT/dp)p = 0) 
from the minima of the pressure along isochoric curves (i.e. (dP/ dT) p = 0). These results are 
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plotted in FigJ7]and are compared with canonical MC estimates. The HMSA is qualitatively 
correct, with errors not exceeding 15 per cent. Nonetheless, once more the inconsistency of 
this LC HMSA is encountered in the results. From a simple thermodynamic analysis it is 
known that the line of maxima of the isothermal compressibility must cross the TMD curve 
at its maximum temperature^. In FigJ7]one can observe that the curve of compressibility 
maxima disappears before reaching the TMD. Moreover, an extrapolation would locate 
the crossing at approximately pa 3 = 0.255 and T* = 0.496, meanwhile the maximum of 
the TMD appears at pa 3 = 0.265 and T* = 0.523. Whereas this can be accepted as 
being qualitatively correct, we find once more that global consistency must be enforced if 
quantitative predictions are to be obtained. 

B. Solid-fluid equilibria 

As explained in Subsection III Bl the complete phase diagram, including the fee solid phase 
which appears at moderate density™, is computed by means of a combination of Wang- 
Landau Monte Carlo simulation (for the determination of LV and LL equilibrium curves), 
NPT MC zero pressure calculations to determine the VS equilibrium curves, standard ther- 
modynamic integration and Gibbs-Duhem integration to evaluate the various fluid-solid 
equilibria. A detailed phase diagram is presented in Figs. [S] and M Note that for densities 
above pa 3 = 0.5 other solid phases are possible. From Jagla's work— (in which the relative 
stability of the zero temperature solid phases is explored) one can infer that at somewhat 
higher densities rhombohedric, cubic and hexagonal (hep) phases could be expected. At still 
higher densities, for which the hard cores play the leading role, again the fee and hep phases 
will be the most stable structures. In any case, from Figure [8] one can already see a fairly 
complex phase diagram resulting from the coexistence of a vapour phase, two liquid phases 
and the fee solid. Inset within Figure M we have enlarged the area of multiple coexistence, 
where one observes the presence of three triple points, namely two neighbouring triple points 
(lower inset) in which one liquid phase (denoted by LI), the vapour and the solid coexist, 
and a third triple point at lower temperature in which the two liquid phases (LI and L2) 
coexist with the fee solid. This latter triple point is also depicted in Figl9] in a detail of the 
P — T phase diagram. 

When the attractive interactions are switched off, even if we have seen that features such 
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FIG. 8: (Colour online) T-p phase diagram of the attractive A2SRP model as obtained from 
computer simulation (circles, LV Wang-Landau MC estimates), (triangles, LL Wang-Landau MC 
estimates), (diamonds, MC NVT thermodynamic integration), (curves, Gibbs-Duhem integration) 
(VS equilibrium line, zero pressure NPT MC). The line connecting the liquid- vapour equilibrium 
points is drawn as a guide for the eye. 

as the curves of compressibility and heat capacity maxima are hardly affected, the phase 
diagram simplifies substantially. Obviously, the LV transition disappears, and one is left with 
a rather peculiar fluid-solid transition (see Fig. [TO]) . First, one observes that when density 
is increased along an isotherm starting from the low density fluid, the systems crystallises 
into an fee phase, which melts upon further compression. At higher density values one could 




FIG. 9: (Colour online) P-T phase diagram of the attractive A2SRP model as obtained from com- 
puter simulation. Points: thermodynamic integration, lines: Gibbs-Duhem integration. Legend as 
in FigEJ 
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FIG. 10: (Colour online) p-T phase diagram of the purely repulsive 2SRP model as obtained 
from computer simulation. Points: thermodynamic integration, lines: Gibbs-Duhem integration. 
Legend as in Fig|8j 

find the rhombohedric, cubic and hexagonal phases predicted by Jagla32. This melting upon 
compression is similar to the behaviour of water. In the 2SRP model this is a purely energy 
driven transition: for certain densities the interparticle distance is necessarily r < d a and 
the repulsive spheres overlap, thus an ordered configuration no longer corresponds to an 
energy minimum. At sufficiently large densities (pa 3 > 0.9), the hard cores will regain 
their controlling role and an entropy driven crystallisation will take place (into either fee 
or hep structures). In the intermediate region (pa 3 > 0.35) the interplay between entropy 
and energy will give rise to a much more complex scenario with different solid phases in 
equilibrium with the fluid phase. 

Finally, from the location of the fluid-solid equilibrium curve, it is now clear that a possible 
LL transition would be preempted by crystallisation. When comparing Figures 151 and [TU1 one 
can see that the effect of the dispersive interactions is to increase the stability of the liquid 
(fluid) phase with respect to the solid. Otherwise, the LL critical temperature happens to 
be fairly close to the the maximum temperature at which the fluid-solid equilibrium takes 
place {T* LL = 0.38 vs T£ s = 0.42). 
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IV. CONCLUSIONS 



In summary, we have presented a detailed study of both the attractive and the purely 
repulsive two scale ramp potential models, using both a self consistent integral equation 
(HMSA) and different Monte Carlo techniques. From this study, it becomes clear that 
the use of a thermodynamically consistent closure is essential if one wishes to capture the 
presence of the LL transition in the attractive model and the density anomaly of the repulsive 
fluid. Moreover, this type of model clearly highlights the shortcomings of local consistency 
criteria, such as those usually implemented in the HMSA. Unfortunately, our calculations 
using a global consistency condition on the pressure calculated from the virial and the 
isothermal compressibility only converge at high temperatures, well away from the LL and 
density anomaly region. This implies, that the closure is missing some essential features 
of the physical behaviour in the region controlled by the soft repulsive interactions. In 
order to bypass this limitation one should explore the use of other functional forms^ or 
use an integro-differential approach of the Self Consistent Ornstein Zernike Approximation 
(SCOZA)£2>^ which has been recently implemented for systems with bounded potentials^. 
Despite these limitations, the theoretical approach yields quantitatively correct estimates 
of the LV transition, the density anomaly and predicts the existence of a LL transition 
(although with a substantial underestimation of the LL critical density). 

By means of extensive Monte Carlo simulations, we have unveiled the rather complicated 
shape of a significant part of the fluid-solid diagram for the A2SRP models, in which three 
triple points have been detected. In the purely repulsive model, one finds that the crystalli- 
sation preempts the liquid-liquid phase separation. This system has, in common with water, 
a solid phase that melts upon compression. In the intermediate region of the phase diagram, 
in which the hard core the the repulsive ramp are competing (i.e. when neither entropy nor 
energy but a subtle combination of both magnitudes leads the system behaviour) one should 
expect even more complex phase diagrams for both attractive and repulsive models. This 
will be the subject of future work. 
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